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ABSt RS i 


A computer model is constructed for the purpose of 
simulating interception of a near earth satellite by a 
ground launched ballistic interceptor with emphasis on 
adaptability of the model to the investigation of interceptor 
performance in terms of interceptor and satellite character- 
istics. The model is keyed to instantaneous satellite 
positions, and provides on call information about minimum 
interceptor velocity requirements, parameters of possible 
interceptor trajectories and relative velocities at inter- 
ception for those trajectories. ‘The model is constructed 
with sufficient generality and flexibility so that it may be 
readily adapted to simulation or computation of other 


intercept problems. 
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The fundamental equations of motion for a body 
acting under the influence of a central force have been 
well known and exhaustively treated since the times of 
Newton and Kepler. These equations, which will not be 
red@rived here, form the basis of the model to be devel- 
oped. It should be noted, however, that this simple 
formulation involves some assumptions that do not com- 
pletely correspond to the physical environment of a near 
earth satellite. These assumptions and their effects on 


the model must be considered. 


BASIC ASSUMPTIONS 

1. The earth and the satellite are considered as point 
masses with the mass of the satellite negligible compared 
to the mass of the earth. 

2. The earth's gravitation is the only force acting on 
the satellite. 


3. The coordinate system is fixed in space. 


DISCUSSION OF ASSUME AAehNS 


If the earth were perfectly spherical with its mass 
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distributed symmetrically about its center, aH its mass 
could be considered to be concentrated at one point. 
However the earth is slightly oblate to a first approx- 
imation and there are slight asymmetries in the distri- 
bution of its mass, particularly in the surface distribution 
of continental and ocean masses. ‘This oblateness is 
responsible for the major anomalies in its gravitational 
field, which affect the behavior of near earth satellites 
in two important ways [3] . First, it causes the orbital 
plane to rotate about the earth's axis in the direction 
opposite to the satellite's motion. This perturbation is 


known as regression of the nodes and its rate is given by 


aS 


SS 
= 9,97 _ ipimeanina cos i degrees/day Cia) 


where R is the radius of the earth, a is the semi-major 
diameter of the orbit, e is the eccentricity of the orbit 
and i is the inclination of the orbital plane to the plane 
of the earth's equator. These terms and symbols are 
discussed in detail in Chapter II. The notation through- 
out this paper will be in accordance with the table of 
symbols unless otherwise specified. Second, the earth's 


oblateness causes the perigee position of the orbit to 








rotate about the principal focus of the orbit with the 
rate being given by 


dW= 4.98 (R_ we weaaa (5 eeesind) degrees/day, (1.2) 


dt a 

The atmosphere offers increasing resistance to the 
motion of a satellite as its altitude decreases. ‘The 
primary effect of atmospheric drag is to reduce the major 
diameter of the orbit, mostly at the expense of apogee 
distance. This makes the orbit more nearly circular with 
the passage of time while reducing the satellite's period 
of revolution and increasing its velocity in accordance with 
Kepler's third law. 

In addition to these major effects, satellite orbits 
are perturbed to a much smaller degree by solar wind, 
radiation pressure and the gravitational attraction of 
other celestial bodies, especially the moon. 

For most applications of the model to be developed, 
the errors mentioned above may be disregarded. The 
model is geared to inspection of instantaneous satellite 
position and velocity vectors, and orbital parameters may 
be considered to be no older than one period of revolution. 


However, the model should not be applied to orbits having 








altitudes of less than 100 miles because of the increased 
effects of atmospheric drag, nor to orbits with apogee 
distances of several earth radii because of perturbations 
due to the moon's gravitational influence. Neither should 
the model be used to predict the position of the satellite 
for more than one or two orbital periods because of the 
cumulative errors indicated by equations (1.1) and (1.2). 
These restrictions are not severe, for in fact most 
objects of interest in the interception problem, such as 
reconnaissance and weapon carrying satellites will orbit at 
altitudes between 100 and 1000 miles due to the require= 
ments of their missions. Substitution of these values 
into equations (1.1) and (1.2) reveal that even in the 
worst cases, the perturbations amount to only a fraction 


of a degree per revolution. 








CHAP Teresi 


ORBITAL MECHANICS 


Kepler's brilliant formulation of the laws of planetary 
motion early in the seventeenth century were provided 
with a sound theoretical basis by Newton's laws of motion 
and were shown to apply generally to bodies acting under 
the influence of a central force. These laws can be 
stated for an earth satellite as follows: 

1. The orbit of a satellite is an ellipse with the 
center of the earth at one of its foci. 

ee The radius vector of a satellite describes equal 
areas in equal times. 

3. ‘The square of the period of a satellite is 
proportional to the cube of its mean distance from the 


center of the earth. 


THE BLLIPTICAL ORBIT 

The size and shape of the elliptical orbit is deter- 
mined completely by the parameters a and e. (See Fig. 
2-1). These quantities can be expressed in terms of 
other orbital dimensions when a and e are not given 


directly. For example, the ellipse is determined if only 
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Figure 2-1. The Elliptical Orbit 
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the minimum and maximum altitudes of the satellite are 

known. ‘The perigee distance p is given by the minimum 
altitude hin plus the radius of the earth R. Similarly, 
the apogee distance q is given by ae + R. The semi- 
major diameter a is then (ptq)/2 and the eccentricity e, 


the semie-minor diameter b and semi-e-latus rectum L are 


given respectively by 


a= —— ae) 

eee a) loc (2-2) 
4 

ee . (2-3) 


from the geometric properties of the ellipse. The angular 
displacement rs) of the satellite from the perigee is 
measured in the direction of satellite motion and is 
commonly called the true anomaly. 

It remains to determine the position of the satellite 
in its orbit as a function of time. In polar coordinates, 
this position is given by 


r = a(lee”) (2-h) 
1te cos 


Unfortunately 9 is not a simple function of time, 


and the satellite position is more easily expressed in 








terms of the eccentric anomaly E, whose geometric 
relationship to the satellite position is illustrated in Fig. 
Ree 

From Newton's laws, the radial acceleration 
(9.80665 m/sec”) of a body in circular orbit of radius R 
is equal to vo /R. Solving for Vo yields a velocity of 
7,904.455 meters/second. 
Then the period in seconds is given by 

Po = 2 = 5104.6536 seconds, (2-5) 
By Kepler's third law, the period of any other near earth 
satellite is then 

a 


B.\2/2 ‘ 
P= Po = seconds, (2-6) 


and the mean angular motion p is given by 

p= att radians/second. (2-7) 
Now define the mean anomaly as follows, 

M =p (t-T) (2=8 ) 
where t is time in seconds and T is the time the satellite 
passes through perigee. From this, Kepler's equation 

M = E -e sin E (2=9) 
can be solved for E. (See Chapter III). From Fig. 2-2, 


it is apparent that 


r= a (l] «e cos EF) (2-10) 
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Figure 2-2. Position vector in terms of E 
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from which, together with equation (2-)) 


cos EB -e 


con® = eS 
1 -e cos & 
(2-11) 
sinO = <Vi1 -e* sin & 


Equations (2-11) make it possible to write the cartesian 


coordinates of S in Fig. 2-2 as 


x' = r cos9@ = a (cos E -e) 
(2=12) 
y' =r sin 9 =aVi1 -e* sin E, 


GEOCENTRIC BQUATORIAL COORDINATES 

The reference system to which the orbital parameters 
are referred is known as the geocentric equatorial 
coordinate system. It is a conventional right-handed 
cartesian system with the origin at the center of the 
earth. (See Fig. 2-3). The x-axis is directed toward 
a fixed point in space where the sun's apparent path 
through the heavens crosses the earth's equatorial plane 
in a northerly direction. This point is variously known 
as the first point of Aries or the position of the vernal 
equinox. It is not truly fixed in space as stated in the 


basic assumptions, but its movement due to nutation of 
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Figure 2-3. Projection of Orbit on Unit 
Geocentric Sphere 
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the earth's spin axis is only a small fraction of a degree 
per year and can be ignored in this problem. The y-axis 

is 90° from the x-axis, measured eastward along the 
equator, and the z-axis is in the direction of the north 
pole. Corresponding unit coordinate vectors are designated 
T, j and k. 

The orientation of the orbit plane in the reference 
system is given by the parameters 1), W, andi. The 
longitude of the ascending node is the angle measured 
eastward along the equator from the x-axis to the point 
where the satellite crosses the equator plane in a northerly 
direction. The argument of the perigee Wis the angle 
measured in the orbit plane in the direction of motion 
from the ascending node to the perigee position. The 
inclination iis the angular separation of the equator and 
orbit planes. This angle is less than 90° if the satellite 
motion is eastward (or direct). If the motion is west- 
ward (retrograde) the angle is the supplement of the 
angular separation of the planes. 

The change from x', y' coordinates in the orbital 
plane to x,y,z coordinates can be accomplished by the 


linear transformation 
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= q ! 
xX P, x' + Q, y 


yr Rx + Q, y' (2-13) 
= q 
Z Ey x' + Q, y 


where the P's and Q's are the direction cosines of the x! 
and y' axes relative to the x-axis [2] . These six 


quantities are found to be 


P, = cos i) cos - sinW sin SL cos i 
Q, = -sinW cos - cos W sin SL cos i 
Py =cosW sinjJ. +sin © cos cos i 
» (21) 
Qy, = -sinW sinfk +cos |W cost cos i 
FZ = sin WJ sin i 
Q, = cos Ww sin i. 


The eccentric anomaly may be introduced directly 


into equations (2-13). If we define 


=—aVi 2 Ow 
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A B 
B..=avVi-e Qy (2-15) 
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then 


a 
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— (cos E -e) + B, sin E 


“< 
ll 


Ay (cos E -e) + By sin EB ; (2-16) 
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(cos E je) + B. sin 
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The first time derivatives of these position 
coordinates yield the components of the satellite velocity 


vector in terms of E. They are found to be 


dx _ 4a : 

== ae (-A, sin E + B, cos E) 

ay . 2 (2 , = 
Te <E (-A, sin E + B, cos E) (2-17) 
dz . we ; 

a <E (-A, sin E + B, cos E) 


CONSHP EVA TION OF ENERGY MOMENTUM hie 
AR BAL, VEEB@Cit ay. 


Kepler's second law and the laws of conservation of 
energy and angular momentum reveal relationships which 
will be found useful. The famous Vis Viva integral, the 


energy equation for an elliptical orbit, states that 


v = GM (2 - = meters”/second”~ (2-18) 


where G is a gravitational constant and M is the mass of 
the earth. In the MKS system of measurement, the 
product GM is equal to 3.98 x 107+ meters”/seconds“ and 
if nautical miles are used as units of distance, GM is 

3 23 : : 
equal to 6.267 x 10° n.m.”~/seconds~. It is noticed 


that for any orbit with a given major axis, v is a function 


Omiy Si re or equivalantly,. of ‘altitude. 


1) 


Angular momentum is conserved in the two body 
problem, and its vector representation is given by the 


constant vector 


+i cones. (2-19) 
dt 


Multiplication by r- yields the equation of the orbital plane, 


r-H =O. (2220) 


From Kepler's second law, the areal velocity of a body in 


an elliptical orbit is constant 2| . kit can be shown that 
a = i, (2-21) 
ait a 


Then by taking the time integral of both sides, the area 


swept out in time t is given by 
Af& = Ht + constant. (2222) 


In particular, if the area swept out between two points ry 


and r, is known, than the time in the orbit between these 


a 


two points is 


Sn Ga Sea 
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CHARA rewiit 


SATERLITE ORD Cert e ttre 


SELECTION OF ORBITAL EBLEMENTS 

A satellite orbit may be determined from a wide 
variety of given parameters. In the present model, 
those initial conditions and parameters are set which the 
analyst may consider to be of particular importance. Of 
first interest are the maximum and minimum altitudes, 
from which the size and shape of the orbit are determined. 
Specifying the inclination immediately determines the 
northern and southern limits of latitude over which the 
satellite may pass and the direction of motion. Setting 
the argument of the perigee fixes the perigee point in 
the desired hemisphere, north or south. Finally, the 
time at perigee and longitude of the ascending node are 
chosen or computed as indicated below and the orbit is 
computed from the equations of Chapter I. 
POSSI TIONRWVEC TORS IN TERMeasOR TERRES liAg 
LATITUDE AND LONGITUDE 

The analyst may wish to know the components of a 


position vector corresponding to a position on earth 


re 








designated by latitude and longitude at a given time. For 
example, he may wish to specify an initial geographic 
position for the perigee of an orbit or determine the 
position vector of the interceptor launch site. Conversely, 
he may wish to know the geographic location of a satellite 
position when its position is given as a vector. 

The Greenwich Hour Angle of the first point of Aries 
(GHA'YT ) is listed in the Nautical Almanac for the begin- 
ning of each day of the year [1] . GHAT is defined as 
the central angle measured westward along the equator 
from the Greenwich meridian to the vernal equinox. 
Designate the longitude of a position by A with eastern 
longitudes being positive, and the latitude by § with 
northern latitudes being positive. Then the Sidereal Hour 
Angle (SHA) between the vernal equinox and the position's 
meridian is equal to GHAT +A »- However SHA increases 
at the rate of the earth's rotation, which in the chosen 
coordinate system is equal to one revolution per sidereal 
day or .0000729 radians/second. If the value of GHAY, 
is taken as that given for oe January 1 of the year, and 
t is the number of seconds after the beginning of the year, 


SHA = GHAY, +A + 0.000729t radians (3-1) 
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and 


r cos SHA cos 5 


> — 
7 = fF sin SHAwcos 6): (3-2) 
Z=r sin $ 


If the position vector r(x,y,z) is given the process 
may be reversed although care must be taken in assigning 
signs to the computed angles. 

z 
= aii (3-3) 
i 2 
§ = (3-1) 
SHA = ain” | hey 33) (3-5) 


where - ee 


" imaeermamnns 


3 


a with the same sign as z and 

- Ti < SHA =7f with its magnitude and sign determined 
by seeing which quadrant contains the x,y coordinates of 
rs \ may then be found directly from equation (3-1). 

If, as suggested above, the analyst wishes to specify 
an initial perigee point and longitude of ascending node from 
geographical coordinates, the corresponding position vectors 
se and fy, are determined from equations (3-1) and (3-2). 
SL is of course simply equal to SHA and 


_ \ 
G = cos : Ea : >} (3-6) 


gee! at 


\ Ss si 


with the same sign as z It must be noted that these 


two points fix the orbital plane so that the inclination may 
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no longer be chosen arbitrarily. However, it is easily 
computed, for the vector 
pM ses OK (3-7) 


is inclined to the z-axis by the angle i if the satellite 


motion is direct and a 5. i if the motion is retrograde. 


vod 
Then 
i (or Ta» oS Ona 





SOLUTIONZOE KE Behari Sy Pepa TIi@eN 

Given the orbital parameters and the satellite's 
position in the orbital plane, the computation of coordinates 
in the geocentric equatorial coordinate system is a straight- 
forward process. However, the position of the satellite 
in its orbital plane as a function of time is found from 
the solution of Kepler's equation which is transcendental 
in EB. ‘This equation is solved by an iterative method of 
succesSive approximations. 

Write equation (2-9) as 

2 omen 2 (3-8) 
and setting E equal to M, solve for EB. Next setting B, 
equal to the computed value of E, repeat the computation. 


The procedure is repeated until the desired accuracy is 


£9 








obtained. It can be shown that the iterated values of 


BE converge rapidly to a single value. 
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CHAP TE Rw 


INTERCEPTOR LRAJHGLG rR. OMe ULATTer 


The first problem to be considered may be stated as 
follows. If an interceptor launch position and a predicted 
intercept point are given, how can one determine trajectory 
parameters that will yield a successful interception. In 
particular, out of all possible interceptor trajectories, 
which one requires the least energy to be supplied to the 
interceptor vehicle. 

After the interceptor has exhausted its propulsive 
energy, it is subject only to the same central force of 
gravitation as is the satellite. It will follow an elliptical 
orbit with a major semi-diameter determined by its altitude 
and velocity in accordance with equation (2-18). Conversely, 
if the semi-diameter of the interceptor orbit is known, 
the energy of the orbit is fixed. Thus the orbit with 
the minimum semi-diameter which passes through the launch 
site and the intercept position is the minimum energy 
interceptor trajectory. 

This reasoning assumes that the interceptor attains 
its maximum velocity immediately after its launch and 


thereafter follows a purely ballistic path with air resistance 


Za 








neglected. Actually there is an initial phase of its flight 
during which it acquires the velocity and attitude 
necessary for the proper ballistic path to the intercept 
position. This phase depends upon the thrust character- 
istics of the propulsion system, the payload weight and 
other factors of interest to the designer and engineer. 
This paper considers only the ballistic requirements which 
the propulsion system must meet, and the interceptor is 
taken to be injected into its orbit with the necessary 
velocity vector at the time of launch. 

Consider the plane determined by the position vector 
of the launch site and the position vector of the intercept 
point as shown in figure -1. 

Let C be the midpoint of line BI and let the 
distance CF= 1/2 GI. (Notation: BI represents the 
magnitude of the directed line segment Bl.) Note that 
by construction the distance OG+GI+tIF = OBt+BE. Now 
define e, as the ratio CI/CF and p, as the quantity 
CF(e,* -i1). Also let Bee = CE (oa Construct the 


hyperbola HH' from the equation 


1 - eh cosQ@ 


Re 











Et 


Fal 


center of the earth Ey 
launch site 
intercept point 


radius of the earth 


satellite altitude 


Figure h-1. Locus of Secondary Foci of 
Interceptor Trajectories 
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with the origin at I. B is the focus of a mirror hyper- 
bola and from the geometric properties of the hyperbola, 
the distances from I and B to any point on the curve 

(E' for example) differ by a constant, in this case equal 
to GI. ‘Thus the distance OB+BF' is equal to the 
distance OI+IF". But recall that the sum of the distances 
from any point on an ellipse to the two foci is a constant, 
equal to the major diameter 2a. Therefore any allipse 
with one focus at the center of the earth and the other 
lying on the hyperbola HH!' will pass through the launch 
site and the point I. Clearly the ellipse with the smallest 
value of a is the one where F'=F", thus determining the 
minimum energy orbit. 

The orbital elements of the minimum energy inter- 
ceptor trajectory may now be found from the given 
position vectors OB and OI. BI=OI - OB. The distance 
GI is the altitude h of the satellite at intercept and is 


equal to OI = Re 


== _ 1 {BIth — 

BF = 5 Goes i (2) 
OF = OB + BE (43) 
a = OB + BEF (md) 


ah 
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e = ° = 
= (4-5) 

The vector product OK = OB x OI is normal to the plane 

of the orbit, and its component will be positive if the 


motion of the interceptor is direct. ‘Then 


da] Cael 





= cose ; =— 
Sz (4-6) 
The vector product ON = k x OK points in the direction 
of the ascending node, so 
ON i 
(ger cor ||P = 
aN (4-7) 


Perigee distance is equal to a - OF, so the perigee position 


vector 


P = (OF -a) OF (4-8) 


and the final spatial element 


—1 | Che. .@ 


(J = cos TON)(OE)|" (4-9) 


The period of the orbit and the mean motion p of 
the interceptor are found from equations (2-6) and (2-7). 
The sine of the true anomaly of the interceptor at the 


launch site is given by 


a5 





sinQ = IO = OB | 


(OP) (OB) egal 
and from equation (2-11), sin E and EF are determined. 
Equation (2-9) gives the value of M and finally 

oo 2a (4-11) 
ya 


from which the epoch of the orbit is determined. A 
similar computation will give the time from perigee to 
intercept from which the time of flight of the interceptor 
is known. 

The above derivation determines the orbital parameters 
for the minimum energy orbit but the same procedure may 
be followed whenever the location of the secondary focus 
is given. ‘The orbital plane is determined by the location 
of the launch site and the intercept point, and any secondary 
focus position in that plane can be computed by equation 
(4-1). A transformation of coordinates will yield the 
position vector of F'' from which the desired parameters 
may be computed by the methods above. 

The velocity vector can now be computed by equation 
(2-17) for any point in the orbit. Of course, the points 


of greatest interest are the point of interception and the 
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launch site. At the intercept point, the difference 
between the satellite and interceptor velocity vectors is 
the closing velocity. ‘The velocity vector at the launch 
site is the sum of the velocities imparted by the propulsion 
system and the velocity vector of the launch site due to 
the earth's rotation. At lower latitudes, this latter 
velocity approaches 1000 mph and an interceptor with a 
given propulsion system will have a substantially longer 
maximum range when launched to the east than when 
launched toward the west. ‘The velocity vector of the 


launch site is seen from equation (3-2) to be 


BES A = cin GleyN ens 1) 
ae 
dy 


aS uR cos SHA cos & (4-12) 


dz 
dt 


This quantity must be subtracted from the orbital velocity 
vector at that point if the velocity to be supplied by the 


propulsion system alone is to be determined. 


Figure 4-2 shows the minimum energy orbit for the 
conditions displayed in Figure 4-1 together with another 
interception orbit having a secondary focus F' on the 


hyperbola HH!', 
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Figure =<. The Minimum Energy Trajectory 
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CHART 2 aay, 


COMPUTER PROGRAMMING @Cehsi DPR te] 


The satellite interception problem, with its complex 
and constantly changing relationships of the many para- 
meters and variables involved lends itself very well to 
simulation by a modern high speed computer. Huge amounts 
of information can be computed, stored and recalled with 
great speed and accuracy. Such a computer model has 
been constructed for execution on a Control Data Corp- 
oration computer model 1604. (See Appendix A). 

Substantial economies can be effected by writing 
computer programs in machine language in terms of 
execution time and memory requirements. This comes at 
the expense of using a format which is not easily adapted 
for use by other computers and is not readily modified by 
other programmers. An alternative, which is the approach 
used in this model,igs to write the computer program ina 
commonly used compiler language, a variation of FORTRAN 
known as FORTRAN-63 |u| - The model follows the 
mathematical approach outlined in this paper, and by 
noting the liberal number of explanatory comments in the 


program, an analyst should be able to adapt the program 


ag 








for his own purposes on another computer employing a 
FORTRAN compiler language. The remainder of this 
chapter is devoted to inspection of computer and compiler 
language considerations which affect the construction of 
the model. 

Inasmuch as the satellite orbit is fixed in the 
chosen coordinate system, much computational time can 
be saved by computing position and velocity vectors for 
one orbital period and storing the information for recall 
as needed. If sufficient storage space is available in the 
computer memory, all quantities of interest such as 
position vector components, velocity vector components, 
true and eacentric anomalies with their trigonometric 
functions, distance from the center of the earth, etc., 
may be stored for later use, and this is recommended 
where possible. However, the period of a circular orbit 
at 1000 miles altitude is nearly 7500 seconds, and storage 
of all quantities of interest in sufficiently small increments 
of time imposes a considerable requirement for memory 
space. 

In this particular model, only the position vector 


components can be stored for each second of time in 


So 








normal format without exhausting the memory storage 
capacity. By going outside the FORTRAN-63 language 
to a compatible assembly program CODAPFP-1, it is 
possible to store the eccentric anomaly for each second 
by storing two quantities in a single storage address [5] ° 
The CODAP-1 subroutine BYTES5 which accomplishes this 
is not included in the FORTRAN-63 listing. 

A small bank of temporary variable names (C1 through 
C9) in a COMMON statement saves some storage by using 
the same names wherever an intermediate, non-permanent 
quantity is computed. The same statement contains 
conversion quantities that are needed in several sub- 
routines. Additional information may be stored on magnetic 
tape or equivalant devices but this is profitable only if the 
time required for access to the information is less than 
the time required to compute it. 

Liberal use is made of subroutines and functions in 
the interest of flexibility, for the model is constructed 
so that more relationships than are explicitly specified may 
be considered with only minor modifications. For example, 
the model can easily consider a terminal guidance aspect 


of the problem or serve as the basis of a probabalistic 


Sid. 





approach to interception. Communication between the 
subroutines is facilitated by the use of COMMON state-= 
ments. 

Clarity of the program for other analysts has been 
enhanced by consistent use of suffixes and prefixes in 
variable names and the use of mnemonic variable names 


throughout. 
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APPENDIX A 
FORTRAN-63 COMPUTER SIMULATION PROGRAM 

The program listed here has been extensively tested 
for a wide variety of orbital elements and the generated 
orbits show close agreement with the orbits of actual 
satellites. Storage is provided for satellite positions in one 
second increments for orbital altitudes to 1000 miles. 
Sample output statements have been included to suggest 
ways in which the user might call for computed values. 

The computational procedure follows the methods 
outlined in the thesis with extensive use of subroutines and 
functions for those calculations that are used repeatedly 
or in other sections of the program. This enhances the 
flexibility of the program for the analyst who may want 
the calculated quantities for supplementary considerations 


im the interception problem. 
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